Thermal transport through a mesoscopic weak link 
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We calculate the rate of energy flow between two macroscopic bodies, each in thermodynamic 
equilibrium at a different temperature, and joined by a weak mechanical link. The macroscopic 
solids are assumed to be electrically insulating, so that thermal energy is carried only by phonons. 
To leading order in the strength of the weak link, modeled here by a harmonic spring, the thermal 
current is determined by a product of the local vibrational density-of-states of the two bodies at 
the points of connection. Our general expression for the thermal current can be regarded as a 
thermal analog of the well-known formula for the electrical current through a resistive barrier. It 
is also related to the thermal Landauer formula in the weak-tunneling limit. Implications for heat 
transport experiments on dielectric quantum point-contacts are discussed. 
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I. INTRODUCTION 



Mesoscopic phonon systems are relatively unexplored 
compared with their electronic counterparts. An excep- 
tion is the recent work on thermal conductance quantiza- 
tion in freely suspended one-dimensional dielectric wires, 
where the thermal conductance was found to be 7rfc|T/6?i 
per transmitted vibrational mode Jl|,^|. This behavior 
parallels the well known electrical conductance quanti- 
zation in units of e 2 /2irh per (spin-resolved) channel in 
one-dimensional mesoscopic conductors |j|-£>)> Electri- 
cal conductance quantization and many other aspects of 
mesoscopic transport in one-dimensional Fermi liquids, 
as well as edge-state transport in integral quantum Hall 
effect systems, can be understood with the Landauer and 
Landauer-Buttiker formalisms |5|||. 

The conventional Landauer formula describes charge 
transport in mesoscopic conductors in the limit where 
there exists one or more propagating channels |Q. An- 
other important transport regime is the weak tunneling 
limit, where the charge conductance is much less than 
e 2 /2irh and, as shown by Schrieffer et al. MM, is deter- 
mined by the density-of-states (DOS) obtained from the 
one-particle Green's function. 

The thermal analog of the weak tunneling limit has 
not been addressed theoretically and is interesting for 
several reasons. First, a microscopic quantum descrip- 
tion of thermal conduction through weak links is crucial 
for understanding energy dissipation in nanostructures 
such as nanoparticles, nanotubes, molecular circuits, and 
nanometer-scale electromechanical systems. As we shall 
demonstrate, the classical theory of thermal conduction, 
based on the heat equation, is entirely inapplicable to 
these systems. Second, thermal conduction through a 
weak link connected to a macroscopic solid turns out to 
be a sensitive local probe of the surface vibrational DOS 
of that solid, suggesting the possibility of a surface mi- 
croscopy based on a scanning thermal probe. 

In this paper we calculate the rate 7th of thermal en- 



ergy flow between two macroscopic bodies, each in ther- 
modynamic equilibrium, and joined by a weak mechani- 
cal link. The weak link may consist of one or more chem- 
ical bonds, or by a narrow "neck" of dielectric material, 
both of which can be accurately modeled by a harmonic 
spring of stiffness K. We obtain a general expression for 
the thermal current that can be regarded as a thermal 
analog of the well-known formula, derived by Schrieffer 
et al. ]§|], for the electrical current through a resistive 
barrier. Our result can also be interpreted as an ap- 
plication of the thermal Landauer formula (||l(J m the 
weak tunneling limit, with the energy-dependent phonon 
transmission probability calculated microscopically. 

Our work is also related to the classic work of Little |Tl] ] 
on the thermal boundary resistance at an interface be- 
tween two dielectrics, a solid-solid analog of the Kapitza 
resistance between solids and superfluid Helium caused 
by phonon scattering at the interface. A tunncling- 
Hamiltonian approach similar to ours has been applied to 
the Kapitza resistance problem by Sheard and Toombs 
]l2[ . In our geometry, however, the thermal resistance 
comes from scattering at the weak link, and the thermal 
current depends on the elastic properties of the link and 
does not vanish if the solids are identical. Heat transport 
in mesoscopic junctions has been studied recently with 
the scattering approach by Cross and Lifshitz |l3| . Ther- 
mal transport through weak links has also been studied 
in conductors, including the two-dimensional electron gas 
|l^ , ^5| and one-dimensional Luttinger liquids jl(J. 

The organization of our paper is as follows: In the 
next section we describe in detail our mesoscopic weak- 
link model, and in Section III we define and calculate the 



local vibrational DOS for the macroscopic solids. A gen- 
eral expression for the thermal current is derived in Sec- 
tion |V|. Some experimental implications are discussed in 
Section [v|, where we calculate the thermal conductance 



through a nanometer-scale junction in Si. Section VI 



con- 



tains a discussion of the differences between electron and 
phonon tunneling, and also of the role of phonon phase 
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coherence in this work. 



II. MESOSCOPIC WEAK LINK 

The model we consider is as follows: Two macroscopic 
solids, L and R, are held at fixed temperatures Tl and 
Tr. The two bodies are assumed to be electrically insu- 
lating, so that thermal energy is carried only by phonons. 
The Hamiltonian of the isolated solids is (we set h = 1) 

H Q = H L + H R , (1) 

where 

Hi = ^w/„a| n a/„, / = L, R. (2) 

n 

The a) j and a n i are phonon creation and annihilation 
operators for the left and right sides, satisfying 

[a n ii a n >i>] = <W &iv (3) 

and 

[a n i, »n'/'] = [a„ 7 , a„/j/] = 0. (4) 

The vibrational modes of the isolated bodies are labeled 
by n and have energies u>i n . Our analysis is valid for any 
LOi n . The mesoscopic weak-link model is illustrated in 
Fig. | 
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FIG. 1. Weak link model. Two macroscopic dielectrics, 
at temperatures Tl and Tr, are joined by a harmonic 
spring of stiffness K. 

The two macroscopic solids are connected by a weak 
mechanical link, which we model by a harmonic spring 
with stiffness K, 

5H = \K(ut~u^) 2 . (5) 



Here uf is the normal component of the displacement 
field at the surface of body / at the point of connection 
to the weak link, with the surface normal taken to be 
in the z direction. The surface displacements can be ex- 
panded in a basis of phonon annihilation and creation 
operators as 

uj = ^ ( h inain + h* In a\ n ) , / = L, R (6) 

n 

where the hi n are model-dependent complex coefficients. 
The hj n appropriate for the stress-free planar surface of a 
semi-infinite isotropic elastic continuum are given below. 

As discussed above, we arc interested in systems where 
the mechanical interaction between the two bodies is ac- 
tually caused by one or more atomic bonds, or by a nar- 
row "neck" of dielectric material. Our harmonic spring 
model correctly accounts for the longitudinal (normal to 
the surface) elastic forces between the solids, but neglects 
any transverse or torsional interaction. Although trans- 
verse and torsional forces can be included by a straight- 
forward generalization of our model, they are often much 
smaller than the longitudinal coupling. 

The macroscopic bodies act as thermal reservoirs, and 
are taken to be ideal thermal conductors. In particular, 
they are assumed to be harmonic [see Eq. (Q)]. Therefore, 
the thermal resistance we calculate is caused entirely by 
the scattering of phonons at the junction between the 
reservoirs and the weak link, and by the finite transmis- 
sion probability through the link. 

III. LOCAL VIBRATIONAL DOS 

In what follows we will need the local vibrational DOS 
(or, more precisely, local spectral density) of the bulk 
solids, evaluated at the point of contact with the weak 
link. These can be obtained from the retarded surface- 
displacement correlation functions 

iM*) = «?(o)]> (7) 

for the isolated macroscopic bodies L and R. Using (||) 
leads to 

Di(t) = -2 9(t) l^"| 2 wn(u In t). (8) 

a 

The local DOS Ni(u>) is then defined in terms of the 
Fourier transform of (Q), 

= -A LnDHw). (9) 

Then we have 

Nj(u) = \hln\ 2 [S(0J - UJln) - 5{lu + u Jn )]. (10) 

n 

In many cases of interest the local spectral density is 
an algebraic function of energy at low energies, 
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JV/(w) cx uj a , 



(11) 



where a is a constant. For example, a — 1 at the planar 
surface of a semi-infinite isotropic elastic continuum (see 
below). 



IV. THERMAL CURRENT 

We now calculate the heat flow between the two bod- 
ies joined by the weak link. The complete system is de- 
scribed by the Hamiltonian 

H = H + SH. (12) 

We define a thermal current operator / t h according to 

I th = d t H R = i[H,H R }. (13) 

The expectation value of 2th is the energy per unit time 
flowing from the left to the right body. 
Writing the interaction (||) as 



SH = \K ^2(A Ln - A Rn )(A Ln , - A Rn 



where 



Aj n = hi n a I„ + h* In a\ 



(14) 



(15) 



we find that the thermal current operator then takes the 
form 



hh — — ^ w Rn { A Rn > — j4l«' , h Rn a Rn — h Rn a Rn } , 

nn' 

(16) 

where {• , •} is an anticommutator. 

The equation of motion for the density matrix in the 
interaction representation is 



d tP (t) = i[p(t),SH(t)], 



where 



Hint /~\^ — iHnt 



0(t) = e %Hot Oe 



(17) 



(18) 



From (0) we find that the nonequilibrium thermal cur- 
rent to leading order is 



Jth(t)=t / dt'([SH(t')J th (t)}) ( 



(19) 



Evaluating Eq. (^9|) leads to our principal result (with 
factors of U reinstated) 

2-kK 2 f°° 

Ith = — ^— J deeN h (e)N R (e) [n L (c) - n R (e)]. (20) 
Here rii,(e) and n R (e) are Bose distribution functions 



n(e) = l/(e £ / fcBT -l) 



(21) 



with temperatures Tl and T R . The details leading to 
Eq. ( 5g ) are given in Appendix 

Our result ( pfj| ) shows that the thermal current between 
a dielectric held at zero temperature and a second dielec- 
tric at temperature T will be a power-law function of T, 
in striking contrast with nonmesoscopic thermal trans- 
port. For example, assuming a spectral density of the 
form (O) leads at low temperature to 



Jth cx T 2 «+ 2 , 

where T is the temperature of the second body. 
The linear thermal conductance, defined by 



(22) 



G t h = km 



(23) 



is given by 



G 



th 



2irK 2 r°° 



dn(e) 



(24) 



This expression, along with Eq. ( |11| ) , shows that the lin- 
ear thermal conductance between two dielectrics held at 
a common temperature T, varies at low temperature as 
a power-law in T, 



G th cx T 



2a + l 



(25) 



where a is the exponent characterizing the power-law 
spectral density at low energies. 



V. THERMAL CONDUCTANCE OF 
NANOMETER-SCALE SILICON JUNCTION 



In this section we give a simple application of our the- 
ory to a structure consisting of a cylindrical neck of Si 
material connecting two semi-infinite Si crystals. To be 
in the mesoscopic regime we assume the dimensions of the 
weak link to be smaller than the phase-coherence length 
of the phonons. The geometry of the system we consider 
is shown schematically in Fig. |[ 
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FIG. 2. Cylindrical silicon junction of length I and di- 
ameter d. 

To apply our formula (^4j) we need the phonon spec- 
tral density at the surface of Si, and also the effective 
spring constant of the link. The spectral density at en- 
ergies much less than the Debye energy may be obtained 
from elasticity theory. This approach, which requires a 
detailed consideration of the vibrational modes of a semi- 
infinite elastic continuum with a stress-free planar sur- 
face, is carried out in Appendix We show there that 
the spectral density at the surface of Si is 



N(e) = Ce, C « 1.3 x 10 8 cm 2 erg" 2 . 
Then using @ we obtain |L7j 

G th = (87r 5 if 2 C 2 /4/15?i)T 3 . 



(26) 



(27) 



The longitudinal stiffness of the mechanical link, a cylin- 
der of length I and diameter d, is 



where Y is Young's modulus. For Si, jli 
Y w 1.3 X 10 12 dyn cm" 2 , 



(28) 



(29) 



and assuming link dimensions of I = 10 nm and d = 1 nm, 
we obtain 



K 



1.0 x 10 4 erg cm 2 , 



and a thermal conductance of 

G th = (9.5 x 10" u erg s" 1 KT 4 ) T 3 . 



(30) 



(31) 



It is interesting to compare this result with the "clas- 
sical" thermal conductance 



G cl - — K 

G th - u K 



of the cylindrical link, as predicted by the heat equation. 
Here k is the experimentally measured bulk thermal con- 
ductivity, itself a function of temperature, which for Si 
can be parameterized as [H 



10 7 erg s 1 cm 1 K 



0.16 + 1.5xl0- 3 T+ 1.6 xl0" 6 T 2 ' 



(33) 



is accurate down 
can use the low- 



wit h temperature in K. Eq. 
to about 100K, below which 
temperature data of Ref . [^0| . 

In Fig. H we plot our result (|3l]) along with the classi- 
cal thermal conductance ([32]) of the same Si link. As is 
evident, these are dramatically different at low tempera- 
tures. The large difference occurs because, as discussed 
m Section [v| the origin of the thermal resistance in the 
two formulas ( |3"l"| ) and ( |32| ) are different. Although it is 
tempting to conclude that they begin to agree at higher 
temperature, this would be incorrect, because our theory 
breaks down at higher temperature and Eq. (|^) is not 
valid up to the temperature where the curves meet. 



^ 0.06 - 




(32) 



300 

T(K) 



FIG. 3. Thermal conductance of the mesoscopic Si link 
shown in Fig. ^, as a function of temperature. The solid 
line follows from Eq. (pSl|), and the dashed line from the 
corresponding classical result given in Eq. (^). The thin 
dotted line is the universal thermal conductance 7rfc|T/6?i 
of a single propagating channel. 

There are four reasons why our analysis becomes in- 
valid as the temperature is increased: The first is that we 
have assumed the weak link to be of mesoscopic dimen- 
sions. As the temperature increases, anharmonic inter- 
action will eventually make the phonon phase-coherence 
length smaller than the size of the link. While an esti- 
mate of the phase-coherence length is beyond the scope 
of this work, the experiment of Schwab et al. [Q suggests 
that in Si it is at least 1 nm at 1 K. The second is that our 
estimate of the spectral density is only valid for tempera- 
tures much less than the Debye temperature of Si, about 
625 K. The third reason is that the leading-order pertur- 
bation theory we have used breaks down when G t h ap- 
proaches the thermal conductance irk^T /6h correspond- 
ing to one propagating channel, shown as a thin dotted 
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line in Fig. [3| And the fourth reason is that we have 
neglected any electronic contribution to Gth, which is 
correct only when k R T is much less than the Si band 
gap. Taking all of these factors into consideration sug- 
gests that Eq. ( |3T| ) is probably not quantitatively correct 
beyond about 10 K. 

VI. DISCUSSION 

In this paper we have studied the thermal analog of the 
weak tunneling limit of charge conduction which might 
be regarded as phonon "tunneling" — and find many sim- 
ilarities to electron tunneling. There are a few important 
differences, however. 

Electron tunneling, as it is usually defined, involves the 
passage of an electron through a classically forbidden re- 
gion. In the thermal case, a phonon of energy e is never 
in a region that does support a mode at that energy. H] 
For example, the harmonic spring employed in Eq. (W) 
can support a single propagating phonon channel, but 
phonons incident on the weak link are mostly reflected 
back into the macroscopic dielectrics. Whereas the tun- 
neling rate of an electron (at a fixed energy e) through a 
forbidden region of thickness / varies exponentially with 



I in the weak-tunneling limit, the thickness dependence 
in the phonon tunneling case is different. In the example 
discussed in Section [v| the thermal conductance varies 
with length I of the bridge as l~ 2 , because the effective 
spring constant (^) of the bridge becomes softer with 
increasing I. 

We have demonstrated that the classical theory of ther- 
mal conduction, based on the heat equation and on the 
concept of a local thermal conductivity, is entirely in- 
applicable to mesoscopic dielectrics. In a mesoscopic di- 
electric, thermal resistance is caused by elastic scattering 
of phonons, whereas in an infinite, disorder-free crystal 
it is caused by inelastic scattering due to anharmonic- 
i ty. In the example of Section M, the quantum result 
( |27| ) is determined by the mechanical properties of the 
bridge material, thr oug h the elastic modulus Y , whereas 
the classical result (|32j ) is determined by the bridge ma- 
terial's bulk thermal conductivity k. 
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APPENDIX A: GENERAL FORMULA FOR THE THERMAL CURRENT 

Evaluating (|l^) we find 

A Lm ,(t') - A Rm ,(t') , (A Ln ,{t) - A Rn ,(t) 



Ith(t) = ^- j dt> W R« Re (l A Un(f) - A Rm (t'), 



x (h K n a Rn (t) - h* Rn a Rn (t)^j 



(Al) 



and, after further simplification, 



Ith{t) = ^T^2I dt' u> Rn Re (( {^w(*Ml„'(0} + {A Rn ,(t),A Rn ,(t')} 



Jo 



h,R n a Rn (t) - h Rn a Rn (t), A Rn (t') 



+ \ h Rn a Rn (t) — h Rn a Rn (t), A Rn (t') \ ( [4„'W^w(t')] + [A R n'(t),A Rn ,(t')] 



(A2) 

where we have used the fact that the commutators are c-numbers. The required thermal expectation values are 

({A In (t), A In (t')}) Q = 2 \h In \ 2 [1 + 2nr(u Jn )] cosw 7n (* - t'), (A3) 



([A In {t),A In (t')]) Q = -2i\h In \ 2 sinuj In {t-t'), 
( [h Rn «Rn (*) - h Rn a f Rn (t), A Rn (i')]> =2| h Rn | 2 cos uj Rn {t-t'), 



(A4) 



(A5) 



and 



({h Rn a Rn (t) - h Rn a Rn (t), A Rn (t')}) Q = -2i \h Rn \ 2 [l + 2n R (uj Rn )] smui Rn (t - t'). 



(A6) 



These lead to 



I th (t) = 2K 2 Y] u Rn \h Rn \ 2 I dt'(\h Lnl \ 2 [l + 2n L (Lu Ln ,)] cosuj Ln ,(t - t') cos co Rn (t - 1') 

+ \h Rn >\ 2 [l + 2n R (u Rn -)] cosuj Rn (t - t') cosui Rn ,(t - t 1 ) 

- \h hn ,\ 2 [l + 2n R (uj Rn )]shiuj Ln '{t-t') smuj Rn (t-t') 

- \h Rn >\ 2 [l + 2n R (uj Rn )] smu Rn (t - t') sinu Rn >{t - t')). (A7) 

Here n^(e) and n R (e) are Bose distribution functions [see Eq. (|2l|)1 with temperatures Xl and T R . Next we make a 
change of variables t' — > t — t' , take the t — *■ oo limit, and include a convergence factor to regularize the long-time 
behavior of the resulting integrals. Finally, using the identities 

dt cosut cosuj't e~ c * = - [S(u - w') + 5(lu + w')] (A8) 



and 



dt s'mujt sino/f e~ c * = - [6(oj - J) - S(uj + u')] , (A9) 



where C is a positive infinitesimal, and reinstating factors 
of h, leads to Eq. (|2C V 

In Eqs. ( po|) and (2^) we have introduced an energy- 
dependent DOS, 

Nj(e) ee ^ \h In \ 2 [S{e ~ hiu In ) - 5(e + kuj In )} (A10) 



which has dimensions of (length) 2 /energy. In a homoge- 
neous elastic continuum of mass density p and volume V, 
N(e) is equal to h 2 /2pe times the thermodynamic DOS 
per volume, V' 1 £) n S(e - e„). 



APPENDIX B: SURFACE DOS OF SILICON 

In this appendix we calculate the local phonon DOS at 
the stress-free planar surface of a semi-infinite isotropic 
elastic continuum, following closely the work of Ezawa 
p2| , and use this to estimate the DOS at the surface of Si. 
The substrate is assumed to occupy the space z > 0. The 
vibrational modes are labeled by n = (to, K, c), where m 
is a branch index taking the values SH, ±, 0, and R, K 
is a two-dimensional wave vector in the xy plane, and 
c ee w/|K| is a parameter (continuous for all branches 
except to = R) with dimensions of velocity. In contrast 
to Ref. we shall use periodic boundary conditions in 
the x and y directions, over a square of area A. 

In our analysis we will approximate Si as an isotropic 
clastic continuum with longitudinal and transverse sound 
velocities 



p = 2.3 g cm 



(B2) 



v\ = 8.5 x 10 5 cm s 1 , 
v t = 5.9 x 10 5 cm s _1 , 



(Bl) 



It will be convenient to treat the Rayleigh branch 
(to = R) separately, and then consider the branches with 
continuous c. In the Rayleigh case the displacement field 
is expanded as [E3J 



K 



2pc R |K| 



[ark fRK + a RK Sk] . (B3) 



where the vibrational cigenfunctions f R K (r) have dimen- 
sions of L~2 and satisfy 

d 3 r • fRK' = <W (B4) 

Here c R — £ vt , where ^ is the root between and 1 of 

f - 8<f + 8(3 - 2v 2 )^ 2 - 16(1 - v 2 ) = 0, (B5) 
and where 

v = v t /v\ (B6) 

is the ratio of transverse and longitudinal bulk sound ve- 
locities. For Si, v = 0.69 and £ = 0.88; hence 



c R = 5.2 x 10 5 cm s 1 . 



(B7) 



The z component of the vibrational eigenfunction at 
the point r = on the surface is 



/r(0) 



2 7 V|K| 



(7 -77) (7-7^ + 27772)^1 



I + 77 2 



and mass density 



(B8) 
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where 



7 = y/i - (c R / Vl ) 2 and V = y/i - (cR/u t ) 2 . (B9) 

We find that the Rayleigh branch contributes to the local 
DOS (H) an amount (for positive u>) 



R branch : 



N(uj) 



gi w 0.42. (BIO) 



Note that <?! generally depends on v, the value quoted in 
( |B10| ) corresponding to Si. 

Next we consider the branches with continuous c. Here 



u = > / dc 



2HKI 



[a m Kc fmKc + fl mKc ClKc] ) 



(Bll) 



where the vibrational eigenfunctions have dimensions of 
L~^c~2 and satisfy 



Kc ' fm'K'c' 



vW<5(c-c'). (B12) 



The range T of the c integration in ( pll ) is [«t,oo] for 
m = SH, [i>i,oo] for m = ±, and [ut,^i] for the m = 
branch. The contribution to the local DOS (for u> > 0) 
from these branches is given by 

N W = T~ E / dc I/™k c (0)| 2 - c|K|). (B13) 

The SH modes are polarized in the xy plane and there- 
fore do not contribute to The ± modes have surface 
amplitude 



/±(0) 



|K| 

AttcA 



±y/a(l + A±iB) 



where 



B 



-7J3 (1 - AtW) . 

= (J3 2 - l) 2 - 4a(3 
~ (/3 2 - l) 2 + 4a/3' 

4VS^(/3 2 -l) 



(/? 2 



4a/?' 



(B14) 



(B15) 



(B16) 



a = y/(c/vi) 2 - 1, and fi = \/{c/v t ) 2 - 1, (B17) 

are all real functions of c. The m = ± branches together 
contribute an amount 



± branches 



N(u) 



Air 2 pvf ' 



flaw 1.0. (B18) 



The value for (72, obtained by doing the integration over 
c in Eq. (B13) numerically, is valid only for the value of 
v corresponding to Si. 

The Til — branch has amplitude 



/o(0) 



IKI 



2ttc/3A 



[- 7 D + i(l -£?)], 



where 



D 



4/3(/3 2 - l) 3 - 16i 7 /? 2 (/3 2 - 1) 



(/3 2 



1) 



16 7 2 /3 2 



and 



E : 



(£ 2 - l) 4 - 16 7 2 /3 2 - 8z 7 /3(/3 2 - l) z 



(/3 2 



1 



16 7 2 /? 2 



(B19) 



(B20) 



(B21) 



are complex-valued functions of c. This leads to a con- 
tribution 



branch 



N{uj) 



87T 2 ptf ' 



53 « 0.59. (B22) 



As before, 53 is obtained numerically and assumes a value 
of f va lid f or Si. Com bining the three contributions 
(plO|), (p!8|), and (p22|), yields |1 



7V(w) 



TlLJ 
47T 2 p 



R 



•1 



93 

2vf 



(B23) 



Using Eq. (|B23|) we obtain the estimate (EQ 
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